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ABSTRACT 

We report preliminary results from a series of numerical simulations of the 
reduced magnetohydrodynamic equations, used to describe the dynamics of 
magnetic loops in active regions of the solar corona. A stationary velocity field is 
applied at the photospheric boundaries to imitate the driving action of granule 
motions. 

A turbulent stationary regime is reached, characterized by a broadband 
power spectrum E\. ~ A;~ 3 / 2 and heating rate levels compatible with the 
energy requirements of active region loops. A dimensional analysis of the 
equations indicates that their solutions are determined by two dimensionless 
parameters: the Reynolds number and the ratio between the Alfven time and 
the photospheric turnover time. From a series of simulations for different 
values of this ratio, we determine how the heating rate scales with the physical 
parameters of the problem, which might be useful for an observational test of 
this model. 

Subject headings: MHD - Sun: corona - turbulence 
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1. Introduction 

Solar coronal loops are likely to be heated by ohmic dissipation of field-aligned electric 
currents. These currents are driven by photospheric motions which twist and shear the 
magnetic fieldlines at the loop footpoints. The generation of small scales in the spatial 
distribution of these currents has been observed and studied in recent simulations (|Mikic| 
et al. 19891 , [Longcope fc Sudan 1994|) . Small scale currents are required by most coronal 



heating theories, since they produce a major enhancement of the energy dissipation rate. 
The development of magnetohydrodynamic (MHD) turbulence in coronal loops has recently 
been proposed ( |G6mez fc Ferro Font an 1992| , [tieyvaerts fc Priest 1992 . Einaudi et al. 1996| , 



Hendrix fc Van Hoven 1996| , |Dmitruk fc Gomez 1997| ) as a likely mechanism to provide such 



enhancement, since the energy being pumped by footpoint motions is naturally transferred 
to small scale structures by the associated energy cascade. 

In the present paper we numerically test this scenario, describing the dynamics of a 
coronal loop through the reduced MHD (RMHD) approximation. In §2 we briefly describe 
the RMHD approximation and perform a dimensional analysis of the equations, which 
leads to an interesting prediction for the functional dependence of the heating rate with 
the physical parameters of the problem. The details of the numerical simulations are 
summarized in §3 and an overview of the results is given in §4. A quantitative scaling law 
for the heating rate is derived in §5, and the relevant results of this paper are listed in §6. 

2. Dimensional analysis of the RMHD equations 

When a coronal loop with an initially uniform magnetic field B = B z is driven by 
photospheric motions at its footpoints, a rather complex dynamical evolution sets in, which 
is described by velocity and magnetic fields. For incompressible and elongated loops, i.e. 
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loops of length L and transverse section (2tt1 ) x (2tt1 ) such that /q ^ L, the velocity 
(v = v(x, y, z, t)) and magnetic B = B z + b(x, y, z, t)) fields can be written in terms of a 
stream function ip = ip(x, y, z, t) (i.e. v = V ± x (ipz,)) and a vector potential a = a(x, y, z, t) 
(i.e. b = V_l x (az)), respectively. The evolution of the fields a and ip is determined by the 
RMHD equations QStrauss 1976|) : 

d t a = v A d g iJ> + [ii>,a]+r}V 2 ± a (1) 

d t w = v A d z j + [i/j, w] - [a, j] + uV 2 ± w (2) 



where v A = Bo/y/Airp is the Alfven speed, v is the kinematic viscosity and r\ is the plasma 
resistivity. The quantities w = —V 2 ± ip and j = — V^a are the z-components of vorticity and 
electric current density, respectively. The non-linear terms are standard Poisson brackets, 
i.e. [u,v] = d x ud y v — d y ud x v. We assume periodicity for the lateral boundary conditions, 
and specify the velocity fields at the photospheric boundaries, 

^ = ) = 0, ^(z = L) = V(x,y) (3) 

where \l/(x, y) is the stream function which describes stationary and incompressible footpoint 
motions. The strength of this external velocity field is therefore proportional to a typical 
photospheric velocity V p . 

To transform Eqs (|I])-(^|) into their dimensionless form, we choose Iq and L as the units 
for transverse and longitudinal distances. Since the dimensions of all physical quantities 
involved in these equations can be expressed as combinations of length and time, let us 
choose t A = L/v A as the time unit. The dimensionless RMHD equations are: 

d t a = d z ij) + [tp, a] + - V' a (4) 

d t w = d z j + [ip, w] - [a, j] + — V \w (5) 
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i 2 i 2 

where S = -j- and R = -®- are respectively the magnetic and kinetic Reynolds numbers. 

Hereafter, we will consider the case S = R, and thus the (common) Reynolds number 

will be the only dimensionless parameter explicitly present in Eqs (§)-© . However, the 

boundary condition (Eqn (H)) brings an extra dimensionless parameter into play, which 

is the photospheric velocity V p divided by our velocity units, i.e. lo/t A . Since V p = lo/t p 

(t p : photospheric turnover time), this velocity ratio is equivalent to the ratio between the 

Alfven time t A and the photospheric turnover time t p . Therefore, from purely dimensional 

considerations, we can derive the following important result: for any physical quantity, 

its dimensionless version Q should be an arbitrary function of the only two dimensionless 

parameters of the problem, i.e. 

Q = J 7 (Qi,Q 2 ), Qi = f f, Q2 = S (6) 

Op 

For instance, let us consider the important case of the heating rate per unit mass, i.e. 
e/p. Its dimensionless version is, 

Q = - n t f = H t f,S) (7) 

P L 1 p 

One of Kolmogorov's hypothesis in his theory for stationary turbulent regimes at very 
large Reynolds numbers ([Kolmogorov 1941]) , states that the dissipation rate is independent 
of the Reynolds number (see also [Frisch 1996|) . In the next section we show that externally 



driven coronal loops eventually reach a stationary turbulent regime. Therefore, hereafter we 
assume that the dependence of the dissipation rate e with the Reynolds number S in Eqn 
(0) can be neglected. For moderated values of the Reynolds numbers as the ones considered 
here, we observed only a mild dependence of the dissipation rate with S (see also |Hendrix 
et al. 1996]) , in accordance with Kolmogorov's assumption. Therefore, 



.i - 



tV K t. 
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We performed a sequence of numerical simulations for different values of the ratio I 1 
to determine the function T . 



3. Numerical simulations 

For the numerical simulations of Eqs (§)-(§|) , ij) and a are expanded in Fourier modes 
in each (x,y) plane (0 < x,y < 2tt and < z < 1). The equations for the coefficients 
ipk(z,t) and a,k(z,t) are time-evolved using a semi-implicit scheme. Linear terms are treated 
in a fully implicit fashion, while nonlinear terms are evolved using a predictor-corrector 
scheme. Also, nonlinear terms are evaluated following a 2/3 fully dealiased (see |Canuto 



|et al. 1988| ) pseudo-spectral technique (see also Pmitruk, Gomez fc DeLuca 1998Q . To 



compute ^-derivatives we use a standard method of finite differences in a staggered regular 
grid (see for instance [Strauss 1976|) . 



We model the photospheric boundary motion in Eqn @ as ^k = ^o = constant 
inside the ring 3 < Hq < 4, and \&k = elsewhere. This choice is intended to simulate a 
stationary and isotropic pattern of photospheric granular motions of diameters between 
27r/o/4 and 27r/o/3. We chose a narrowband and non-random forcing to make sure that 
the broadband energy spectra and the signatures of intermittency that we obtained (see 
below) are exclusively determined by the nonlinear nature of the MHD equations. The 
normalization factor ty Q is proportional to t A /t p , as mentioned above. 

4. Development of MHD turbulence 

We performed a sequence of numerical simulations of Eqs (^)-(^|) with 192 x 192 x 32 
gridpoints, S = 1500 and different values of the ratio t A /t p in the range [0.01,0.15]. The 
typical behavior of the heating rate as a function of time, is shown in Figure 1 for the 
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particular case of t A /t p = 0.064. After an initial transient, the heating rate is seen to 
approach a stationary level. This level is fully consistent with the heating requirements for 
coronal active regions of 10 7 erg cm~ 2 seg" 1 flWithbroe fc Noyes 1977|) . The intermittent 



behavior of this time series, which is typical of turbulent systems, is also ubiquitous in 
all our simulations. Once the stationary regime is reached in each of the simulations, we 
determine the mean value and rms of the departure of the series from the mean. 

A stronger indication of the presence of a turbulent regime, is the development of 
a broadband energy power spectrum, which behaves like E^ ~ k~ 3 ^ 2 for both two and 
three dimensional MHD turbulence (Kraichnan spectrum). To compute the energy power 
spectrum, we performed a single simulation with better spatial resolution (384 x 384 x 32), to 
allow for a more extended inertial range. The parameter for this simulation is t A /t p = 0.064 
and S = 5000. The spectrum shown in Fig. 2 was taken at t = 5t p , i.e. well in the 
stationary regime. This higher spatial resolution is still insufficient for the formation of a 
broad inertial range to determine the slope of these spectra with confidence. However, the 
observed slope is consistent with a Kraichnan spectrum, i.e. E& — k~ 3 ^ 2 . Also, note that 
the kinetic and magnetic power spectra reach equipartition at large wavenumbers, even 
though the total kinetic energy is much smaller than the total magnetic energy. 

Another important consequence of the energy power spectrum displayed in Fig. 2, is the 
fact that smaller scales dissipate more energy than the larger scales, since e& oc k 2 E^ ~ k x l 2 . 
Three dimensional simulations of loops driven by random footpoint motions have shown 
that energy dissipation preferentially occurs in current sheets which form exponentially fast 
( pVlikic et al. 1989J , |Longcope fc Sudan 1994J ). Within the framework of MHD turbulence, 
two dimensional simulations also show the ubiquitous presence of current sheets evolving in 
a rather dynamic fashion ( [Matthaeus fc Lamkin 1986J , [Biskamp fc Welter 1989] , pmitruk" 



Gomez fc DeLuca 1998 ). The presence of elongated current sheets as the most common 



dissipative structures is also apparent in these simulations (also |Hcndrix fc Van Hoven| 
|1996| ). Figure 3 shows the spatial distribution of the electric current density j(x,y,z,t) 
for our 384 x 384 x 32 simulation (t A /t p = 0.064 and S = 5000). Intense currents going 
upward (downward) are shown in white (black) at various planes z = constant, and for 
t = 5t p . Although the reconnection in these turbulent regimes is expected to be fast ( Priest 



fc Forbes 1992| , [Hendrix fc Van Hoven 1996| , |Milano et al. 19§9"D , the Reynolds numbers 



attained by current simulations cannot provide a definite answer to this question. 

5. Determination of the heating rate scaling exponent 

One of the main goals of this paper is to derive the functional dependence of the 
heating rate with the ratio between the two relevant timescales (see Eqn (|8D). To this end, 
we performed various simulations for different values of the parameter t A /t p . Figure 4 shows 
that this functional dependence can be adequately fit by a power law, i.e. 

e = f^r, 8 = 1.51 ±0.04 (9) 

Each diamond corresponds to the mean value of the heating rate for a particular 
simulation in its stationary regime. The error bars correspond to the rms value of the 
departure from the mean. The slope in this plot corresponds to the parameter s quoted in 
Eqn (H), which was derived following a minimum squares procedure. 

This result is fully consistent with the prediction arising from a two dimensional MHD 
model ( |Dmitruk fc Gomez 1997|) , s 2D = §, which is assumed to simulate the dynamics of a 



generic transverse slice of a loop. This coincidence between both scalings suggests that two 
dimensional simulations provide an adequate description of the dynamics of coronal loops 
(see also |Hendrix fc Van Hoven 1996J ), provided that the external forcing in Eqs (]l|)-(@) is 
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riven by v A d z ip « -j- 2 - and v A d z j ~ ( [Einaudi et al. 199C 



Dmitruk fc Gomez 19971) 



6. Discussion and conclusions 

In the present paper we report preliminary results from a series of numerical simulations 
of the RMHD equations driven at their footpoints by stationary granule-size motions. The 
important results are summarized as follows: 



(1) After a few photospheric turnover times, the system relaxes to a stationary turbulent 
regime, with dissipation rates consistent with the heating requirements of coronal 
active regions (10 6 — 10 7 erg cm~ 2 s^ 1 ). 

(2) A dimensional analysis of the RMHD equations shows that the dimensionless heating 
rate does only depend on the ratio t A /t p and the Reynolds number. However, for 
turbulent systems at very large Reynolds numbers, the heating rate is likely to be 
independent of the Reynolds number. From a series of numerical simulations, we find 
that the heating rate scales with t A /t p ase~ ttKt 1 ) 3 '' 2 , i-e. 



pi(B V p ) j ! i 



L 



L' 



(10) 



which might be useful for an observational test of this theoretical framework of coronal 
heating. 

(3) The dissipative structures observed in our simulations are current sheets elongated 
along the axis of the loop. The spatial and temporal distribution of these structures is 
rather intermittent, as expected for turbulent regimes. We associate these dissipation 
events taking place inside coronal loops to the nanoflares envisioned by [Parker 1988| . 
as a likely scenario for coronal heating. 
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Fig. 1. — Dimensionless heating rate as a function of time (in units of t A ) for t A /t v = 0.064 
and S = 1500. The thin trace corresponds to viscous heating. 

Fig. 2. — Energy power spectrum for a 384 x 384 x 32 simulation at t — 5t p . The thin trace 
corresponds to kinetic energy. A straight line corresponding to a slope —3/2 is displayed for 
reference. 

Fig. 3. — Halftones of the function j(x,y) at the transverse planes z = 0,0.33,0.66, 1 and 
for t — 5t p . Upflowing currents are white and downflowing currents are black. 

Fig. 4. — Dimensionless heating rate as a function of t A /t p . Each diamond corresponds to a 
different simulation and the full line corresponds to the slope determined from a minimum 
squares fit. 
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